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INTRODUCTION 



I . 

Two important parameters of interest in the field of aero- 
dynamics of airfoils are lift and drag. These can be evalu- 
ated either experimentally or theoretically. The desire for 
computational methods to aid the design process is promoted 
by the need to reduce the number and cost of wind tunnel tests. 

This thesis treats the problem of incompressible, two- 
dimensional steady flow about airfoils or airfoil combinations 
at large angles of attack. Such flows exhibit strong viscous 
flow effects including regions of flow separation. Therefore 
methods are required which can account for these effects. 

Currently there exist two main methods , namely the direct 
computation of viscous flows by means of the Navier-Stokes 
equations or the so-called viscous/inviscid interaction method. 
The former approach is more straightforward but also much more 
expensive and time-consuming. Therefore, the latter approach 
is to be preferred if it can be shown that it produces good 
agreement with the available experimental results. 

It is the purpose of this thesis to contribute to the 
evaluation of the viscous/inviscid interaction method. To 
this end, the viscous/inviscid computer codes developed by 
Cebeci and collaborators at the Douglas Aircraft Company were 
obtained and applied to several airfoils. 

In addition, a separate panel method was formulated and pro- 
grammed in order to obtain the inviscid flow over airfoil 
combinations . 
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The basic equations are formulated in Chapter II. Chapter 
III addresses the problem of inviscid flow calculations using 
the panel method. In Chapter IV the solution of the boundary 
layer equations by means of the Cebeci-Keller box method is 
explained. Finally, Chapter V describes the viscous/inviscid 
interaction problem and presents results of computations. 
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II. BASIC EQUATIONS 



A. INTRODUCTION 

In this chapter, the basic equations of fluid flow are 
derived. We find that the resulting equations are PDE ' s whose 
exact solutions exist only in very few cases. The PDE ' s are 
classified, "parabolic," "hyperbolic," "elliptic" depending 
on the geometry of their domains of dependence and regions of 
influence, and the solution procedures are different according 
to the type of equation. Table 2-1 gives a brief classif ication 
of these equations. 



TABLE 2-1 



CLASSIFICATION OF PDE'S 



1 Elliptic 



Parasol ic 



Hyperbolic 






Physical Upstream 

Meaning Influence 



No Upstream 
Influence 



No Upstream 
Influence 



Example • Laplace 

Equati ons 

• Steady Navier- 
Stokes 



• Thin Shear 
Layer 



• Supersonic 
Flow 



: P is perturbation point 




is domain on which solution 
at P depends 




is region of influence of 
a perturbation at P 
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For example, the Laplace equation is elliptic. Its solution 
would have to be repeated for many iterations so that the up- 
stream influence can be gradually propagated (panel method in 
Chapter III), but the thin shear layer equations are parabolic. 
Their numerical solution can be obtained by marching in the 
downstream direction only (Box method in Chapter IV) . 

B. DERIVATION OF GENERAL EQUATIONS 

The continuity equation and the Navier-Stokes equations 
are basic for an aerodynamic analysis. We start with the basic 
physical concepts and derive the general equations for 2-D, 
unsteady, compressible, viscous flow. 

1 . Continuity Equation 

One of the basic laws of "Newtonian mechanics" states 
that mass can neither be created nor destroyed. Therefore, 
for a fixed control volume (see Figure 2.1), the principle of 
mass conservation can be stated that the net mass flow rate 
into and out of the control volume equals the time rate of 
change of mass within the control volume. 

If the central point ' P' has representative fluid proper- 
ties (velocity, density, pressure, etc.), then properties at 
other locations can be obtained by Taylor series expansions. 
Therefore the x-component of the velocity at the center of the 
positive x-face (right-hand face) is 



u 



+ 



3u ,dx, 
3x 2 



+ 



3x 



u ,dx, 2 1 
2 1 2 ‘ 21 



+ . . . 



( 2 . 1 ) 
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Figure 2.1. Control Volume for 2-D 



As dx goes to zero, all higher order terms will be dropped, 
so that only the first two terms will be considered. Similarly, 
the density at the center of the positive x-face is 



P 



+ 



9 p dx 
3x 2 



( 2 . 2 ) 



Then the mass flow rate out of the positive x-face is 



Mass Flow Rate Out 



(Density) (Velocity) (Area) 






w , 8u dx , 

. . . ) (u + 8x 2 



[pu + |_(pu)3£]dy 



. . . ) dy 

(2.3) 



By the same procedure the mass flow rate into the control volume 
through the negative x-face (left-hand side) is 



16 



Mass Flow Rate In 


= tpu - f^(pu) ^]dy (2.4) 


From Eqs. (2.3) and 


(2.4) , we get the net mass flow rate 


through the control 


volume in the x-direction. 


Net Mass Flow Rate 


= tpu - |^-(pu)^-]dy - tpu + |^(pu)^rldy 




= - |^(pu) dx dy . (2.5) 



In a similar manner, the net mass flow rate in the y-direction 
is 





- 1— (pv) dx dy (2.6) 

dy 


The total mass flow 


rate through the control volume is 


obtained by summing 


Eqs. (2.5) and (2.6). 


Total Net Mass Flow 


Rate = - [f^r(P u ) + f^( pv ) ] dx dy (2.7) 



Next, we consider the time rate of change of mass within 
the control volume. 

Time Rate of Change 3 , , , . 

of Mass = 3t (p dxd y> 

= |#dxdy (2.8) 

o t 
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Now we combine Eqs . (2.7) and (2.8) using the concept of 



conservation of mass. Then 




-[ fx (pu) + f^-(Pv) ]dx dy = |^dxdy 


(2.9) 


Therefore we obtain the general form of the 
equation for two-dimensional flow as 


continuity 


3(pu) 9 ( pv) 9p _ n 

9x 9y 9t 


(2.10) 


For steady or unsteady incompressible flow, 


Eq. (2.10) reduces to 


+ *2L = o 

r\ 1 r\ v • 

9x 9y 

2. Navier-Stokes Equations 


(2.10a) 



Newton's second law of motion, when mass is conserved, 
equivalently states that the rate of change of the momentum of 



a body equals the sum of the forces applied 


to that body, or 


T~» H 

l f = at (mv) 


(2.11) 



In considering a small volume element of fluid, there are two 
types of forces to be considered, namely surface forces which 
are acting on the surface and body forces which are acting on 
the fluid inside the elemental volume, such as gravity (see 
Figure 2.2). 
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Figure 2.2. 



Forces Acting on the Fluid 



Assuming that the stresses are known at point 'P', we get 
expressions for the stresses on the fluid element surfaces 
by Taylor series expansion. 



Net force in 
x-direction due 
to Normal Stresses 



3a 



[a + • 
xx 



xx dx. , r 

t -=-] dy - [a 

dX 2 2 XX 



da , 
xx dx, , 

— — ld y 



3d xx . , 

~xr dx d y 



(2.12) 



Net force in 




x-direction due 


= [a 


to Shear Stresses 


yx 




3a 




yx 




3y 


Therefore the total 


surface 


formed by summation 


of Eas. 



3a , 

_^x dy ]dx . [ 0 

ay 2 yx 




dv 

2 



] dx 



dx dy 



(2.13) 



forces in the x-direction are 
(2.12) and (2.13) . 
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J f (surface) 



(2.14) 



3a 



= [- 



xx 



3a 



yx-i 



3x 



3y 



dxdy 



Let f (body) be defined as the body force per unit mass with 
the following components : 



? (body) = X i + Y j 



Thus, the x-component of the body force is 



f (body) = m X 

= p dxdy • 1 - x 



(2.15) 



Adding Eqs. (2.14) and (2.15) provides the total force in the 
x-direction . 



I F 



f (surface) and f (body) 

X X 



3a 



= [pX + 



xx 

3x 



3a 



+ “^r ] dxdy 



(2.16) 



Now we consider the rate of change of the momentum of the fluid. 

Let us take the x-component only. Then, since the mass is constant 




x 



m 



du 

dt 



= p dxdy (u 



3u 

3x 



+ v 



3u 3u 
3y 3t 



(2.17) 
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because, u = f [x (t ) ,y (t) , t] , and by chain-rule 



du _ 3u dx 3u dy 3u 3t 

dt 3x dt 3y dt 3t 3t 



= u 



3u 

3x 



+ v 




3u 
3 1 



Substitution of Eqs . (2.16) and (2.17) into the x-component 

of Eq. (2.11) gives the final equation of motion. 



pX + 



3o 



xx 

3x 



+ 



3o 



yx 

3y 



r 3u 

pIu 53? 



+ V 



3u 3u. 
3y 3 1 



(2.18) 



Now we want to show the stress in terms of the velocity com- 
ponents. In this thesis we will consider only simple "Newtonian" 
fluids obeying Stokes' law. This means that the 'extra' stress 
(above hydrostatic pressure) is proportional to the rate of 
strain. With the definition, 



Extra Stress = constant x (rate of strain) 



and introducing y = coefficient of viscosity. 



a 



xx 



- p + 2 p <n> 



xy 



,3u , 3v. 

+ 35' 



where the pressure in an incompressible fluid is seen to be 
equal to minus one-third the sum of the three normal-stress 
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components in view of Eq. (2.10a), In two dimensions then, 

a + a = -2P . Eq. (2.18) then becomes, if the body force 
xx yy ^ 2 

is neglected 



3u 

3t 



+ 



u 



u 

3x 



V 



3u 

3y 



t t 3 o' •n 3 o' 

i i xx + 1_ xy 

p 3x p 3x p 3y 



(2.19) 



Substitution then produces, for incompressible flow. 



3u 
3 1 



+ 



u 



u 

3x 



v 



3u 

3y 



1 _3P 
p 3x 



+ v[ 




+ 



,2 
3 Ui 



3y' 



( 2 . 20 ) 



where v = y/p = kinematic viscosity, and similarly for the 
y-component 



3 v 
3 1 



+ 



u 



3v , 3 v 

3x 3y 



if* + v t i^ + 

p 3x 2 



»2 
3 v- 

2 - 



( 2 . 21 ) 



These are the well-known Navier-Stokes equations for two- 
dimensional incompressible viscous flow. 



C. INVISCID FLOW EQUATION 

All real fluid flows are viscous, but inviscid flow can 
be assumed outside of a thin boundary layer and a narrow wake 
behind the body for large Reynolds numbers. This is the reason 
why the inviscid flow equations are important even though they 
represent an ideal case. If the flow far upstream is uniform 
then it is also irrotational . This allows the introduction of 
the velocity potential. 
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1 . Velocity Potential 



The velocity potential iji is a function whose gradient 
represents the fluid velocity. Thus 



34 > 

3x 



u 




( 2 . 22 ) 



where 



<p = <p (x,y , t) 

Therefore, the significance of the velocity potential lies 
in the fact that one equation for <p can be used rather than 
three equations for the velocity components. 

2 . Laplace Equation 

For steady, incompressible flow, the continuity equa- 
tion (2.10) becomes 



3u 3_v 
3x 3y 



0 



(2.23) 



This equation can be written in terms of velocity 
potential <p by substituting Eq. (2.22). Thus 





0 



(2.24) 



This is the well-known Laplace equation (vector form is V (ji = 0) . 
It is a linear equation which makes it possible to apply the 
principle of linear superposition. For instance, 
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2 

If <j>^, <p^, . .., <J> n satisfy V <J> = 0 , then also 

2 

4> = <J>^ + <j >2 + ... + <J) satisfies V <p = 0 . 

Thus the flow resulting from the superposition of incompressi- 
ble, irrotational flows is also an incompressible and irro- 
tational flow. This superposition principle makes it possible 
to build up quite complicated flow from a few simple solutions 
of Laplace's equation. The singularity (or panel) methods 
presented in the next chapter are based on this idea. 

D. THIN SHEAR LAYER EQUATIONS 

High Reynolds number flows over airfoils (and other con- 
figurations) generate a very thin shear layer (boundary layer) 
close to the airfoil surface. If 6 denotes the boundary layer 
thickness and x the distance from the leading edge of a flat 
plate, then it is well-known that 

6 (x) ~ /vx/U or ^ ~ /1/Re 

X X 

where v is the kinematic viscosity and Re = Ux/v. This formula 

X 

shows that 



« 1 if Re » 1 

X X 

Hence the flow outside of the boundary layer can be considered 
to be inviscid, but the effect of viscosity cannot be neglected 
within this layer. Nevertheless, the Navier-S tokes equations 
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for steady incompressible flow can be simplified because of 
the fact that 6 is much smaller than the representative length 
of the airfoil (the chord length) . This can be deduced from 
the Navier-Stokes equations as follows: 



1 

p 


+ 

3 x 


a 2 

, 3 u 
V ( — ~ + 

3 x 


~2 

3 u x 

2 ’ 

3 y 


9 u 

= u 33 F + 


9 u 

v 

9 y 


1 

p 


3 y 


a 2 

, 3 v 
V ( p + 

3 x 


~2 

3 v s 

2 ' 

3 y 


9 v 

= U 3 ^ + 


9 v 

v 

3 y 



u 


is 


now 


replaced 


by 


a 


typical 


value , 


say 


u 

e 


X 


is 


now 


replaced 


by 


a 


typical 


value , 


say 


l; 


y 


is 


now 


replaced 


by 


a 


typical 


value , 


say 


6 . 



Then ^ can be expressed by U e /6; 

can be expressed by U^/Z; 

d P 9 

-r — can be expressed by p\J z /Z 
ax e 



( 2 . 25 ) 



( 2 . 26 ) 



(because P and U e are related by the Bernoulli equations) . 
Therefore the x-component terms of the Navier-Stokes equation 
can be estimated to have the following magnitudes: 



1 

P 



3P 

3x 



a 2 i 2 

, a u . 3 U, 

V ( T + T ) 



3x 



3y 



u 




V 



3u 

3y 





u 



2 

e 



i 
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The magnitude of the term v ^ follows from the application 

° y u 

of the continuity equation = 0, i.e., — ^ ~ = 0 or 

dx dy o 

v ~ ^ and therefore 



v 



9u 

9y 



U U 
r e e 

5 T T 




The two viscous terms are of vastly different magnitude 

because 3^u/9x^ ~ U /Z^ and 3^u/3y^ ~ U^/6 ^ hence 

2 2 2 2 2 2 
3 u/3y >> 3 u/3x and 3 u/3x can therefore be neglected 

2 2 2 2 
compared to 3 u/3y . Finally, the term v 3 u/3y must be of 

the same magnitude as the other terms if the influence of 

viscosity is to be retained. The y-component terms of the 

Navier-Stokes equations are easily estimated to be smaller 

than the x-component terms because 



u 



3v 

3x 



U 



_ 6 _ 

,2 





6 



and hence are smaller by a factor j-. Therefore the two 
equations reduce to 



u 





1 3P 
p 3x 



+ v 




(2.27) 



3P 

3y 



o 



(2.28) 



By adding the continuity equation (Eq. (2.23)) to these rela- 
tions, we get the basic equations to describe laminar flow 
thin shear layers. 
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E . TURBULENT FLOW 



We must deal with the instantaneous properties in the 
turbulent flow. Thus the time-mean value concept is applied: 

u = u + u ' 

where u is the time-mean value, and u' is the fluctuating 
portion. Similarly, 



v = v + v ' 



P = P + P' 



Introducing these relations into Eq. (2.20) 



u 



9u 

3x 



+ 




1 8P ^ , 9 



+ 



3y‘ 



r) u 



9u 1 u 1 _ 9u ' v 1 
9x 9y 



(2.29) 



We can see that pu'u' and pu'v' correspond to a normal stress 
and a shear stress. We call these terms turbulent stresses 
or Reynolds stresses. 

Similar analyses can be done for the y-component equations 
and z-component equations in the three-dimensional case. The 
extra Reynolds stresses can be summarized by the following 
array , 
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a 


a 


a 






u 1 


1 ^ 


u 


'v' 


u 


'w' 


XX 

a 


xy 

a 


XZ 

a „ 


_ 


-p 


u 1 


'v' 


V 


,2 


u 


1 w 1 


yx 

a zx 


yy 

a zy 


yz 

a zz 
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Much of the effort in turbulent flow studies centers on the 
proper modelling of these turbulent terms. 
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III. PANEL METHOD 



A. INTRODUCTION 

The panel method was developed in the 1960's at McDonnell 
Douglas by Smith and Hess as a numerical approach for 2-D 
and 3-D potential flow problems. This chapter presents the 
application of the panel method to 2-D problems about one or 
two bodies. The basic idea consists in representing the flow 
past a body by a distribution of singularities (sources, 
sinks, vortices) on the surface such that the body surface 
becomes a streamline. 

The numerical approach requires some approximations (the 
assumption given in parentheses refers to our approach) . 

A. The surface of the body is replaced by a finite 
number of elements (straight-line-elements) . 

B. The condition of tangential flow is satisfied at a 
finite number of points, the so-called control- 
points (midpoints of elements). 

C. The singularity distribution of each element is 
approximated by some kind of analytical functions 
(singularity strengths are assumed to be constant 
along any one element, but vary from element to 
element) . 

The advantages of the panel method in comparison with other 
procedures are: 

A. The panel method does not include an approximation in 
the physics — thin airfoil theory does. 

B. The panel method can be easily applied to both 2-D 
and 3-D problems — a virtually unsolvable task for 
conformal mapping procedures, which are confined 
to 2-D configurations only. 
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C. The panel technique can be readily extended to flow 
fields past several bodies--a task which causes at 
least some troubles in "classical" mapping techniques. 

The method's versatility has been proven in various extensions, 
e.g., hydrofoils, cascades, nozzles, and even complete air- 
craft. Since its origin the method has been improved by 

features like higher order approximations to both body surface 
and singularity distribution, taking account for compressi- 
bility effects, and inclusion of wake models. Today the panel 
method is probably our most powerful tool in analyzing poten- 
tial flows. 

B. NONLIFTING FLOW PAST A BODY 

The effects of lift (respectively, camber and angle of 
attack) and displacement (resp. thickness) can be studied 
separately because of the linearity of Laplace's equation. 

This section is concerned about displacement flows due to the 
thickness of bodies, a flow which is usually represented by 
sources and sinks. 

We will first draw the reader's attention to a single 
straight-line-element, along which sources of constant strength 
are distributed. This simple case allows us to explain the 
basics of the panel technique. The source strength A is de- 
fined as the volume of fluid discharged per unit area. Since 
fluid is ejected perpendicular to the panel's surface in both 
directions, discharge velocities are half of the source 
strength. The boundary condition for an inclined panel requires 
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that the normal component of the free stream velocity is 
balanced by this discharge velocity (see Figure 3.1) . 

J - - cos 6 (3.1) 




Figure 3.1. Boundary Condition at One Inclined Panel 

This relation between a geometric quantity (3) and the unknown 
source strength establishes tangential flow on the panel 
surface . 

Things, which are obvious for a single panel, become 
slightly more complicated for a structure of panels. Mutual 
interference of source panels, i.e., each panel induces a 
velocity at other panels, causes the complication. While the 
boundary condition of a single panel had been set up by glancing 
at a simple geometry sketch, we now better switch to a syste- 
matic procedure, emphasizing the concepts of velocity poten- 
tial and superposition. We consider a 2-D closed body. 
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approximated by several panels and inclined at an angle to 
the oncoming flow. Our goal is to derive a relation for the 
unknown source strengths from the condition of tangential flow 
at the control points. To this end we will give the velocity 
potentials of a single source, a source panel, and a closed 
body built up by a source panel, in that order. 

Radial streamlines and concentric circular equipotentials 
characterize one of the very basic potential flows, the single 
source flow. Its velocity potential is defined by 



where r is the distance from (x,y) to the source. Arranging 
single sources on a straight line element corresponds in 
terms of the velocity potential to a summation of single poten- 
tials, in the limit of an infinite number of sources to an 
integration over the panel length. Thus the velocity potential 
of a source panel can be written as 



m source panels, representing a body, induce a flow field, 
whose velocity potential at any point (x,y) is given by 



single source <{>(x,y) 




( 3 . 2 ) 



source panel c))(x,y) 



7r- f Jin r dS 
2tt > _ 



( 3 . 3 ) 



♦ (x<y) = l <(i.(x,y) 

3=1 3 



m 




Jin r . dS . 
3 3 



( 3 . 4 ) 
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Figure 3.2. Designations for Calculation 



We call 4> the potential of the flow disturbances due to the dis- 
placement flow. The total velocity pptential of a nonlifting 
flow results from a superposition of this displacement flow 
and a uniform flow, which is inclined at an angle a to the 
x-axis . 



Recalling the definition of the velocity potential (velocity 
equals gradient of potential), the boundary condition of tan- 
gential flow takes the form 



$(x,y) = V cosax + V sin ay + 7 — 1 

0O OQ m *-» 



, a . 

m A . j 



j=l 



£n r . dS . 
3 3 



3S> 

3n 



0 on the surface 



Applying this condition within the framework of the panel 
method provides a system of equations , 
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j = l 



for i = 1 , . . . ,m 



(3.5) 



which establishes zero normal velocity at all control points. 

This linear system can be solved for the unknown source 

strengths after the integrals have been evaluated. 

Concept of influence coefficients 

The influence coefficient, I. is defined as the normal 

iD 
t h 

velocity at the i panel due to a source distribution of 

th 

2iT-strength at the j panel. 
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D i 



(3.6) 



th 

The contribution to the normal velocity at the l panel by 

t h 

the actual source distribution of the j panel is the product 
of A j/2n and the influence coefficient I^j . To compute the 
influence coefficient we must substitute 



r ij = V < x Mi ' X j )2 + <y M . 'V 2 
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in Eq. (3.6) and carry out the differentiation. 
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where 



9x . 
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9n. 

1 



= cos 3. 



9Yi 
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sin 3 • 

l 



x . = x_ + S . cos 0 . 

3 B j 3 3 



y . = y + S. sin 9. 

3 B j 3 3 



t h 

The integration co'vers the length of the j panel. Finally 
we obtain 



[x M - (x B +S .cos0 . )]cos6 i + [y M - (y B +S . sin0 . ) ] sin3 i 
J i i 3 -* j j 3 



1 • • = / ' 1 2 2 2 

13 0 [ * M -(x B +S.COS0.)] +[y -(y B +S.sin9 )] 

ij JJ ij JJ 



ds . 

3 



(3.8) 



Equation (3.8) is expressed in terms of only and, after 
some manipulations, the integral may be written in the form 



b - cs . 

I ±j = / "2 

13 0 S -eS . +f 

3 3 



(3.9) 



where 



b = (X M " X B )C ° S 3 i + (Y M." Y B. )sin 6 i 
*i j 1 3 



c = cos 0jCos 8^ + sin 0^ sin 6^ 
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e = 2 cos0.(x M -x ) + 2 sin6.(y M -y ) 

i j i j 

f = <X M.' X B . )2 + 

1 D ID 

The integration can now be performed analytically. 
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Determination of unknown source strengths A 

Equation (3.5), expressed in terms of Eq. (3.9) and divided 
by V,, takes the form 



m 

tt A ! + J All.. = - cosacosg.- sin as in g. 

1 j=l D ID i i 



(3.11) 



where 



A'. 

D 



A . 

1 

2ttV 



or in the more convenient matrix form 
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36 



The above set of linear equations can be solved for A^ by 
Gauss' elimination method or any other linear equation 
algorithm. 

On-body velocities 

t h 

The velocity at the midpoint of the i panel, Vg can 

l 

be obtained by a spatial derivative of the velocity potential 
in tangential direction. 
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(3.14) 



where J . . = / |g— (Jin . ) dS . denotes the tangential velocity 

13 j i ID D 

at the i th panel due to a source distribution of 2iT-strength 
fch 

at the j panel . 
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The calculation of J\ ^ follows the same procedure as that 



of I . . , so 
13 




J. . = 

13 


0 

j b -cs. 

= / -5 3 dS. (3.15) 
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Positive signs of on-body velocities indicate that velocities 
are oriented in the direction of the surface coordinates, 
while negative signs indicate opposite directions of veloci- 
ties and surface coordinates. The positive direction of 
surface coordinates is defined clockwise. Therefore positive 
values of on-body velocities are to be expected on the upper 
surface, negative values on the lower surface. 

Off-body velocities 

Streamlines can be determined by computing velocities at 
off-body points and using a numerical quadrature to progress 
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from one point to another point on a streamline. The velocity 
components can be expressed according to 



U(x ± ,y i ) = 



8 $ (x i ,y i ) 
3 oT 
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Normalizing the above equations by the free stream velocity 
and abbreviating the integrals simplify the relations to 
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I X . and iY . are again influence coefficients, whose evaluation 
ID ID y 

can be adopted from the already introduced procedure. 
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where 
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C. CIRCULATORY FLOW 

While inviscid 2-D flow theories are not capable of pre- 
dicting drag characteristics, information about lift can be 
provided by them. Creation of lift is closely related to a 
type of flow called circulatory flow. We mentioned already 
that the flow around a lifting airfoil can be decomposed into 
two elementary flows, i.e., displacement flow and circulatory 
flow. Circulation and circulatory flows are the subjects of 
this section. 

The early approaches of airfoil theory emphasized a flow 
model in which the airfoil was represented by an infinitely 
thin vortex sheet only. This so-called thin airfoil theory 
predicts lift quite well, because lift depends primarily on 
circulatory flow. Unfortunately a straightforward extension 
of vortex sheets to "surface singularity" method is impossible. 
Therefore aerodynamicists have proposed a couple of flow models. 
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which allow the implementation of circulatory flows in 
"surface-singularity" methods. Examples are: 

(1) Smith and Hess represent circulatory flows by a com- 
bination of source and vortex distributions. [Ref. 1] 

(2) Martensen prefers vortex distributions only, but states 
the problems in terms of the stream function. [Ref. 2] 

(3) Davenport makes use of linearly varying vortex 
distributions. [Ref. 3] 

Our approach follows the ideas of Smith and Hess. These circu- 
latory flows are composed of a vortex distribution, which is 
constant along all and for all panels, and a source distribution 
of conventional shape. 

We start at the very beginnings of vortex flows. Concentric 
circular streamlines and radial equipotentials characterize the 
flow field of a single vortex. Its velocity potential can be 
written as 



single vortex c})(x,y) 



_T_ 

2ir 



arctan 



y-y 
J J v 

x-x 

V 



(3.20) 



with (x ,y ) as the center of the vortex. A structure of m 
vortex panels induces a flow field, whose velocity potential 
at any point (x,y) is given by 

£ . 

m n i y-y^ 

4* (x , y ) = r l (- j-) f arctan — — 1 dS . (3.21) 

j=l Z1T 0 X j D 

This flow field differs in two important points from the non- 
lifting flow field: 

(1) It violates the condition of tangential flow. 

(2) The unknown singularity strength F cannot be 
determined immediately. 
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The task of determining circulation must be postponed to the 
implementation of the Kutta condition. Temporarily we set 
the vortex strength equal to one. Tangential flow must be 
established by the aid of an additional source distribution. 
Strengths of this additional source distribution must be com- 
puted according to the condition that the normal velocities 
due to the vortex distributions at the control points are 
balanced by the normal velocities due to the additional source 
distribution . 
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Abbreviating the integrals by the above defined influence 
coefficients, we get 
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where 



A (1) 

1 



are the unknown strengths whose effect is intended 
to balance the normal velocities induced by a unit 
vortex distribution. 




is the normal influence coefficient due to a source 
distribution . 




is the normal influence coefficient due to a vortex 
distribution . 




is the tangential influence coefficient due to a 
source distribution. 




is the tangential influence coefficient due to a 
vortex distribution. 



42 



Since influence coefficients of source and vortex distributions 



(v) 



are related by I . . 

ID 

according to 
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= -J\ j , the above equation can be expressed 
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(3.23) 



By solving this system for , we determine the properties 

of circulatory flow of unit strength. 

Calculation of disturbance velocity due to unit circulatory 
flow 

The disturbance velocity, , is composed of two parts, 

one due to the constant vortex distribution, the other due to 
the additional source distribution 
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Making use again of the relation between influence coefficients 

( J. (v) = iff^), we have 
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D. SYNTHESIZING A COMBINED FLOW 

The Kutta condition serves as matching condition for nonlift- 
ing and circulatory flow. These two basic flows must be 
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superimposed such that flows of upper and lower surface merge 
smoothly at the trailing edge. This original version of the 
Kutta condition is usually substituted by the condition of 
zero load (or equal velocities on both upper and lower sur- 
face) at the trailing edge (see Figure 3.3). 




Figure 3.3. Single Airfoil: Superposition of Nonlifting 

and Circulatory Flow, Controlled by the 
Kutta Condition 

Since the panel method does not permit the evaluation of 
velocities at the trailing edge, the Kutta condition is satis- 
fied approximately by requiring that velocities at the control 
points of the rearmost panels have equal magnitude. Therefore 
the rearmost panels should be chosen short so that flow at 
their midpoints will effectively represent that at the trailing 
edge . 

Determination of circulation 

Suppose 1 S ^ and m 1 "^ panels are the closest panels to the trail- 
ing edge on the lower and upper surface, then we can write the 
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Kutta condition as 





V (N) 

m 



+ r v 



(V) 

m 



(3.26) 



where VV denotes the tangential velocity in nonlifting flow. 

Equation (3.26) can be solved for the circulation T. 

Calculation of on-body velocities and of pressure coefficients 

Three parts contribute to the total velocity: free stream, 

disturbance due to displacement flow and disturbance due to 

(N) 

lifting flow. Say V designates the velocity due to the 
nonlifting flow including the free stream component and 
represents the velocity due to a lifting flow of unit circu- 
lation. Then the total tangential velocity at the midpoint of 
th, 

the i panel is given by 



V. = vf N) + Tvf v) (3.27) 



Once the velocity has been computed, the pressure, customarily 

expressed by means of a dimensionless coefficient C , is 

F 

determined by Bernoulli's equation: 



C 




V. , 

1 - (— ) 2 

CO 



(3.28) 



Addendum : More than one body configuration. 

One of the main advantages of the panel method is its 
easy extension to multi-element airfoils. As a matter of 
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fact even flow past an infinite number of bodies can be solved 
by means of the panel method, if these bodies are arranged in 
form of cascades. The minor changes, which are necessary to 
apply the panel method to a finite number of bodies, include: 

(1) The overall scheme must provide a circulatory flow 
for each lifting body. (The number of nonlifting 
flows remains one.) 

(2) Flow past each body with lift is subject to a Kutta 
condition. Accordingly the numbers of equations 
requiring zero load at the trailing edge equals the 
number of circulatory flows, which allows the 
definite determination of each lifting body's 
circulation . 

Figure 3.4 illustrates the superposition of nonlifting- and 
circulatory flows for a two element airfoil. 




Figure 3.4. Two Element Airfoil: Superposition of 

Nonlifting and Circulatory Flows 



46 



E . EXAMPLES 



This section illustrates the capabilities of the program 
"PANEL" which can be applied to 2-D potential flow problems 
past one or two bodies. 

Flow past one circular cylinder 

The source panel technique is applied to the flow past 
a circular cylinder. This case is regarded as nonlifting, 
i.e., the cylinder experiences no force perpendicular to the free 
stream. As sketched in Figure 3.5, the surface of the cylinder 
is approximated by eight panels of equal width. For zero angle 
of attack, Eq . (3.5) reduces to 

A . m A. 3 ( S,n r . . ) 

V cos0. +~+ l f — yr dS . = 0 (3.29) 

00 i 2 . u , 2 tt ' . d n . i 

3 = 1 3 i 

Solving a set of 8 simultaneous algebraic equations, the source 
strengths and the pressure coefficients can be determined. 
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Figure 3.5. Arrangement of Panels on a Circular Cylinder 
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The results are shown in Figure 3.6 where they are com- 



pared with analytical results (C =1 -4sin 6) 

P 
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Figure 3.6. Pressure Coefficient on a Circular Cylinder 
Obtained by Using Eight Source Panels 
(Marked by 0) in Comparison with the Exact 
Solution 



This example demonstrates the power of the panel method. 
However the reader should be aware that only 8 panels are not 
sufficient to describe the geometry in most of the cases. 
Basically the achieved accuracy depends on both the shape of 
the body and the panel configuration (number of panels and 
local widths) . A closer spacing is advisable in regions where 
severe changes of the pressure distribution are expected 
(e.g., leading edge). 
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Flow past a pair of circular cylinders 



Two circular -cylinders are arranged side-by-side in a uni- 
form stream. The surface of each cylinder is replaced by 50 
panels of equal width (see Figure 3.7(a)) . The computed 
velocity distribution on one of the cylinders is shown in 
Figure 3.8(b) . The reader shall pay some attention to a 
comparison between the flow past one and the flow past a pair 
of cylinders. Obviously the maximum velocity is increased 
by the existence of a second cylinder. The closer the two 
cylinders are arranged, the higher the maximum velocity. 

While the stagnation points in a single cylinder flow are 
located at the farthest down and upstream points of the cylinder, 
the disturbance of a second cylinder causes the stagnation 
points to move towards the other cylinder. The streamline 
picture, given in Figure 3.8, should provide a deeper under- 
standing of this kind of flow. 

Flow past two element airfoil 

The main goal of leading and trailing edge devices is to 
obtain a higher lift coefficient. We will investigate the 
effect of a single slotted flap on the pressure distribution 
of the main airfoil. 

The pressure distributions of a single airfoil and of an 
airfoil-flap combination are compared in Figure 3.9. The 
coordinates of both main airfoil (a NACA 4412) and flap are 
listed in Section F (sample input data) . The results indicate 
that lift increases more than 50% by using a slotted, 21.5 
degree deflected flap. 
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0C- 





(b) 

Figure 3.7. (a) Arrangement of Panels on Two Cylinders 

Side by Side 

(b) Calculated Velocity Distribution on One 
of Two Identical Circular Cylinders Whose 
Centers Are One and a Quarter Diameters 
Apart 
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Figure 3.8. Streamlines of Flow Past Two Identical Cylinders 
Whose Centers Are One and a Quarter Diameters 
Apart 
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x/c 

Figure 3.9. Pressure Distributions on a Single and 
a Flapped Airfoil 



Airfoil in ground effect 

Flow past an airfoil in ground effect is another applica- 
tion of our program. 

The boundary condition at the ground requires vanishing 
normal velocity there. We meet indirectly this condition by 
arranging the second, imaginary airfoil such that the ground 
becomes an axis of symmetry of this "new" flow field (see 
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Figure 3.10). Since an axis of symmetry must be impermeable 
to fluid particles, the desired flow is obtained without 
explicitly satisfying the boundary condition at the ground. 

This kind of flow is a challenge to aerodynamicists for several 
reasons. Whenever an airplane takes off and lands, it passes 
a zone where flow is severely affected by the proximity of the 
ground. Wind tunnel experiments must be corrected for wall- 
effects, quite a similar situation with grounds below and above 
the airfoil. And there was a German experimental seaplane 
that makes use of flying very close to the sea level. However, 
our numerical experiments will tell only one part of the story 
because all these flows are highly 3-dimensional. 



Real 




Ground 



Image 



Figure 3.10. Airfoil in Ground Effect 

Let's first question how does the pressure distribution 
change near the ground. Figure 3.11 shows that lift is reduced 
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Figure 3.11. Pressure Distributions on a Single NACA 4412 
and on a NACA 4412 in Ground Effect 
(h/c = 0.2, a = 5°) 
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on the upper surface and increased on the lower surface. In 
the particular case the overall lift gain is about 15% of the 
lift in free air, but we might not always expect a lift gain. 
The actual balance between lift reduction on the upper and 
lift increase on the lower surface depends on both distance 
from the ground and angle of incidence. Figure 3.12 confirms 
that there are cases with less lift than in free air. High 
angles of attack and moderate distances from ground are sus- 
ceptible constellations to lift loss. 

F. I/O — DESCRIPTION AND LISTING OF THE PROGRAM "PANEL" 

This program calculates non-lifting and lifting potential 
flow past one or two bodies. Any 2-dimensional shape and 
any angle of attack, which do not cause flow separation, 
are acceptable. 

Input data 

The data must be arranged in the following order: 

(1) Header card; 

(2) Coordinates of first body cards (variable number of 
cards) ; 

(3) Second body control card; 

(4) Coordinates of second body cards (variable number of 
cards) . 

Items 3 and 4 are used only for the 2-body case. The actual 
instructions are as follows. 

Header card 

1-10: Number of bodies (integer) 

11-20: Number of points of the first body (integer) 

21-30: Angle of attack in degrees (real) 
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Coordinates of first body cards 



The input procedure of body coordinates requires the follow- 
ing sequence: start at the trailing edge, progress on the 

lower surface to the leading edge, return on the upper surface 
to the trailing edge and finish with the trailing edge. The 
trailing edge of closed bodies is input twice, as first and 
last point. However airfoils with finite trailing edge thick- 
ness should be treated as non-closed bodies, i.e., the last 
point input is not the first point repeated. 

1-10: X coordinates of the points defining the body (real) 
11-20: Y coordinates of the points defining the body (real) 
Second body control card 
1-10: Number of points of the second body 

Coordinates of second body card 

The X- and Y-coordinates of the 2 nc ^ body are input in the 
same format as the coordinates of the first body. 

Output 

There are two kinds of solutions, non-lifting and lifting, 
both of which are preceded by the following column header. 

PANEL X Y V C 

P 

where 



PANEL is the number of the panel; 

X and Y are the coordinates of control points (not 
boundary points) ; 

V denotes the relative velocity (V/V ) ; and 
C denotes the pressure coefficient. 
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Sample problem 



This sample illustrates program input and output. The 
data refer to the airfoil-flap example of Section III.E 
(see Figure 3.9). 
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Input 



1 . 



1 1 

1 . 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 
1 . 



35 10 . 

.0 

-.0016 

-.0022 

-.0039 

-.0065 

-.01 

-.014 

-.018 

-.0226 

-.025 

-.0274 

-.0288 

-.0286 

-.0274 

-.0249 

-.0195 

-.0143 

.0 

.024 

. 0339 

.0473 

.0576 

.0659 

. 0789 

.0880 

. 0941 

.0976 

.0980 

.0919 

.0814 

.0669 

.0489 

.0271 

. 0147 

.0 

-.15 

-.14 

-.125 

-.08 

-.05 

-.06 

-.09 

-.115 

-.15 



X . . .1 

2 

95 

90 

80 

70 

60 

50 

40 

30 

25 

20 

15 

10 

075 

05 

025 

0125 

0 

0125 

025 

05 

075 

10 

15 

20 

25 

30 

40 

50 

60 

70 

80 

90 

95 

0 

9 

25 

20 

15 

05 

0 

05 

15 

20 

25 
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Output 



NONLIFTING SOLUTION 



PANEL 


XM 


YM 


V 


CP 


1 


0.97500 


-0.00080 


-3.05132 


-8.31055 


2 


0.92500 


-0.00190 


-1.71788 


-1.95110 


3 


0.85000 


-0.00305 


-1.43382 


-1.05583 


A 


0.75000 


-0.00520 


-1.20962 


-0.46318 


5 


0.65000 


-0 .00825 


-1.11906 


-0.25229 


6 


0.55000 


-0.01200 


-1.07063 


-0.14625 


7 


0 .45000 


-0.01600 


-1.03543 


-0.07210 


8 


0.35000 


-0.02030 


-1.01792 


-0.03616 


9 


0.27500 


-0.02380 


-0.97360 


0.05210 


10 


0.22500 


-0.02620 


-0.98044 


0.03874 


11 


0.17500 


-0.02810 


-0.97498 


0.04942 


12 


0.12500 


-0.02870 


-0.95981 


0 . 07876 


13 


0.08750 


-0.02800 


-0.92944 


0.13613 


14 


0.06250 


-0.02615 


-0.90437 


0.18212 


15 


0.03750 


-0.02220 


-0.81613 


0.33393 


16 


0.01875 


-0.01690 


-0.66716 


0.55489 


17 


0.00625 


-0.00715 


-0.06604 


0.99564 


18 


0.00625 


0 .01200 


1.20372 


-0.44893 


19 


0.01875 


0.02895 


1.51674 


-1.30050 


20 


0.03750 


0 . 04060 


1.46366 


-1.14230 


21 


0.06250 


0.05245 


1.45104 


-1.10550 


22 


0 . 08750 


0 . 06175 


1.44313 


-1.08263 


23 


0.12500 


0.07240 


1.40582 


-0.97632 


24 


0.17500 


0.08345 


1.38729 


-0.92458 


25 


0.22500 


0.09105 


1.36530 


-0.86405 


26 


0.27500 


0.09585 


1.34794 


-0.81694 


27 


0.35000 


0.09780 


1.24546 


-0.55117 


28 


0.45000 


0.09495 


1.15851 


-0.34214 


29 


0.55000 


0.08665 


1.07036 


-0.14566 


30 


0.65000 


0.07415 


0.96331 


0.07203 


31 


0.75000 


0.05790 


0.82090 


0.32612 


32 


0.85000 


0.03800 


0.55160 


0.69573 


33 


0.92500 


0.02090 


0.37152 


0.86197 


34 


0.97500 


0 .00735 


-0.86208 


0.25681 


35 


1.22500 


-0 .14500 


-2.69341 


-6.25447 


36 


1.17500 


-0.13250 


-1.67415 


-1.80278 


37 


1 .10000 


-0.10250 


-0.59324 


0.64806 


38 


1.02500 


-0.06500 


1.34805 


-0.81725 


39 


1.02500 


-0.05500 


3 . 04627 


-8.27976 


40 


1.10000 


-0.07500 


1.49617 


-1.23852 


41 


1.17500 


-0 .10250 


0.51932 


0.73031 


42 


1.22500 


-0.13250 


-1.15598 


-0.33629 



59 



1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 



LIFTING SOLUTION 



XM 


YM 


0.97500 


-0 .00080 


0.92500 


-0.00190 


0.85000 


-0.00305 


0.75000 


-0 .00520 


0 .65000 


-0 .00825 


0.55000 


-0.01200 


0.45000 


-0.01600 


0.35000 


-0 . 02030 


0.27500 


-0.02380 


0 . 22500 


-0.02620 


0.17500 


-0.02810 


0.12500 


-0.02870 


0.08750 


-0.02800 


0.06250 


-0.02615 


0.03750 


-0.02220 


0.01875 


-0.01690 


0.00625 


-0.00715 


0.00625 


0.01200 


0.01875 


0.02895 


0.03750 


0.04060 


0 . 06250 


0.05245 


0.08750 


0.06175 


0.12500 


0.07240 


0.17500 


0.08345 


0.22500 


0 . 09105 


0.27500 


0.09585 


0.35000 


0 . 09780 


0.45000 


0 . 09495 


0.55000 


0.08665 


0.65000 


0.07415 


0.75000 


0.05790 


0.85000 


0.03800 


0.92500 


0 . 02090 


0.97500 


0.00735 


1.22500 


-0.14500 


1.17500 


-0.13250 


1 .10000 


-0.10250 


1.02500 


-0.06500 


1.02500 


-0 . 05500 


1.10000 


-0 . 07500 


1 .17500 


-0.10250 


1.22500 


-0.13250 



V 


CP 


-1.09525 


-0.19958 


-0.80975 


0.34430 


-0.71356 


0.49083 


-0.65850 


0.56638 


-0.63122 


0.60157 


-0.60472 


0.63432 


-0.56660 


0.67897 


-0.52361 


0.72583 


-0.47182 


0.77738 


-0.43188 


0.81348 


-0.35514 


0.87387 


-0.22992 


0.94714 


-0.07829 


0.99387 


0.09794 


0.99041 


0.47846 


0.77108 


1.09475 


-0.19848 


2.42874 


-4.89880 


3.46087 


-10.97764 


3.13743 


-8.84344 


2.71609 


-6.37716 


2.48044 


-5.15257 


2.34667 


-4.50686 


2.20878 


-3.87871 


2.08988 


-3.36759 


2.00707 


-3.02833 


1.94596 


-2.78674 


1.81808 


-2.30540 


1.69455 


-1.87149 


1.60093 


-1.56299 


1.51477 


-1.29451 


1.43499 


-1.05920 


1.33681 


-0.78706 


1.28943 


-0.66263 


1.09525 


-0.19956 


-0.86324 


0.25482 


-0.79007 


0.37579 


-0.43515 


0.81064 


0.41219 


0.83010 


1.91951 


-2.68453 


1.49954 


-1.24862 


1.23353 


-0.52159 


0.86324 


0.25482 
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Program listing 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 



c 

c 

c 

c 

c 

c 

c 

c 



THIS PROGRAM CALCULATES 2-D POTENTIAL 
AT ANY ANGLE OF ATTACK. 



FLOW PAST 1 OR 2-BODIES 



WRITTEN BY 
DATE 



CAPT LEE, HEE WOO 
NOV. 28 1985 



NOTE i MAXIMUM NUMBER OF PANELS = 200 



C 

C 

C 

C 

C 

C 

C 



C 

C 

C 



xccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DIMENSION Z( 20 0,200) , XB( 200 ) , YBC 200 ) , BE( 200 ) ,THC200) ,V(200) 
DIMENSION WKAREAC 65000), XMC 200 ) , YM( 200),VVV(200),VC(200) 

DIMENSION CP (200) , VV( 200 ) , Y( 200 , 200 ) , AC 200 ) , VI ( 200) , SI ( 200 ) 
DIMENSION VC1C200), VC2 ( 200 ) , V2( 200 ) , VTC200) 



READ INPUT DATA (FOR FIRST BODY) — 

READ( 4 , 1 ) NB, NN, AN 

I FORMAT (2110, F10. 5) 

AN = ANX3. 141592/180. 

DO 10 I = 1 , NN 

READ(4,11)XB(I),YB(I) 

II FORMAT ( 2F1 0.5) 

10 CONTINUE 



C 

C 

C 



CALCULATE MID-POINTS OF PANELS AND ANGLES THETA 



N = NN-1 
DO 12 1=1, N 
K = 1+1 

XMC I ) = (XB(I)+XB(K) )/2. 

YM(I) = (YB(I)+YB(K) )/2. 

THH = (YB(K)-YB(I))/(XB(K)-XB(I)) 

TH( I ) = ATAN(THH) 

IF(XB(K) .LT.XB(I)) TH( I ) =TH( I ) +3 . 141592 
12 CONTINUE 
C 

C CALCULATE PANEL LENGTHS 

C 



DO 13 I = 1 , N 
K = 1 + 1 

A( I ) = SQRT((XB(K)-XB(I))XX2+(YB(K)-YB(I))X*2) 

13 CONTINUE 

IFCNB.NE.l) GO TO 600 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c c 

C NON-LIFTING PART (1-BODY) C 

c c 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO 14 I = 1 , N 
C 

C — - CALCULATE ANGLE BETA AND NORMAL COMPONENT OF FREE STREAM 
C VELOCITY — 

C 



BE(I)= TH(I)+3. 141592/2. 

VCI) = -(C0S(BE(I))XC0S(AN)+SIN(BE(I))XSINCAN)) 

DO 15 J = 1,N 
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C CALCULATE INFLUENCE COEFFICIENTS OF NORMAL VELOCITY 

C 

IFCI.EQ.J) GO TO 15 

B = CXMCI)-XBCJ))xCOSCBECI))+CYMCI)-YBC J ) ) XSINC BEC I ) ) 

C = COSCTHCJ))XCOSCBECI))+SINCTHCJ))XSINCBECI)) 

E = 2.XC0SCTHCJ))X(XMCI)-XBCJ))+2.xSINCTHCJ) ) xC YMC I)-YB(J)) 
F = CXMCI)-XBCJ))xx2+CYMCI)-YBCJ))xx2 
H = SQRT(4.XF-EX*2) 

ZCI,J) = TEGCB,C,E,F,H,ACJ)) 

15 CONTINUE 
ZCI,I)=3. 141592 

14 CONTINUE 
C 

C SOLVE SET OF LINEAR EQUATIONS FOR SOURCE STRENGTHS 

C 

IDGT = 0 

CALL LEQT2F ( Z, 1 , N, 200 , V, IDGT, WKAREA , I ER) 

C 

C CALCULATE INFLUENCE COEFFICIENTS OF TANGENTIAL VELOCITY 

C 

DO 16 I = 1,N 

DO 17 J = 1 , N 

IFCI.EQ.J) GO TO 17 

B = CXMCI)-XBCJ))XCOSCTHCI))+CYMCI)-YBCJ))xSINCTHCI)) 

C = C0S(TH(J))XC0S(TH(I))+SIN(TH(J))XSIN(THCI)) 

E = 2.XC0S(TH(J) )X(XM(I)-XB( J) )+2.XSIN(THC J) )X(YM(I)-YB( J) ) 
F = CXMCI)-XBCJ))xx2+CYMCI)-YBCJ))xx2 
H = SQRTC4.XF-EXX2) 

YCI,J) = TEGCB,C,E,F,H,AC J)) 

17 CONTINUE 

Y(I, I) =0.0 

16 CONTINUE 
C 

C CALCULATE TOTAL VELOCITY AND CP AT MIDPOINTS OF EACH PANEL 

C 

WRITEC 8 ,95) 

95 FORMATC/// , 25X, 'NONLIFTING SOLUTION*,// 

X10X, ’PANEL * ,5X, *XM* ,8X, *YM* , 11X, *V* ,10X, *CP* ) 

DO 18 I = 1 , N 
S = 0 . 

DO 19 J = 1 , N 

S = S+VC J ) XY( I , J ) 

19 CONTINUE 

VV(I) = COS(TH(I))xCOS(AN)+SIN(TH(I))xSIN(AN)+S 
CP(I) = 1.-VV(I)XX2 
WRITEC 8, 93) I , XMC I ) , YMC I ) , VVC I ) , CPC I ) 

93 FORMAT C10X,I3,3X,2F10.5,3X,2F10.5) 

18 CONTINUE 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C LIFTING PART C 1-BODY) 

c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C CALCULATE SOURCE STRENGTHS DUE TO CIRCULATORY FLOW OF UNIT 

C STRENGTH 

C 

DO 20 I = 1 , N 
S = 0. 

DO 21 J = 1 , N 

S = S+Y C I , J ) 

21 CONTINUE 

V1CI)=S 

2G CONTINUE 
IDGT = 0 

CALL LEQT2F CZ, 1 , N, 200 , VI , IDGT, WKAREA, I ER ) 
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non 



ono o o o non odd non no no 



— CALCULATE DISTURBANCE VELOCITY DUE TO CIRCULATORY FLOW OF UNIT 

STRENGTH 



DO 22 I = 1 , N 
S = 0. 

DO 23 J = 1,N 

S = S+Y (I,J)XV1(J)+Z(I,J) 

23 CONTINUE 

VCCI)=S 

22 CONTINUE 

CALCULATE VORTEX STRENGTHS BY KUTTA CONDITION — 

SI = -( VV( 1 )+VV( N) )/(VC(l)+VC(N)) 

CALCULATE TOTAL VELOCITY AND CP AT MIDPOINTS OF EACH PANEL — 

WRITEC 8 , 96 ) 

96 F0RMAT(///,27X, 'LIFTING SOLUTION',// 

XI OX, 'PANEL ',5X, 'XM',8X, 'YM',11X, 'V • , 10X, 'CP ' ) 

DO 24 I =1 , N 

VVV(I) = VV(I)+SI*VC(I) 

CP(I) = l.-VVV(I)x*2 

WRITEC8 , 93)1, XM(I),YM(I), VVVCI ) , CP( I ) 

24 CONTINUE 

25 FORMAT ( 3F1 0.5) 

GO TO 700 

READ INPUT DATA (FOR SECOND BODY) 

600 READ( 4,31) MM 

31 FORMAT (110) 

DO 32 I = NN+1 , NN+MM 

READ(4,11)XB(I),YB(I) 

32 CONTINUE 

M = NN+MM-2 

CALCULATE MID-POINTS OF PANELS AND ANGLES THETA — 

DO 33 I=N+1,M 
K = 1 + 1 

XM( I ) = (XB(K)+XB(K+1) )/2. 

YM( I ) = ( YB ( K) +YB (K+l) )/2. 

THH = (YB(K+1)-YB(K))/(XB(K+1)-XB(K)) 

TH( I ) = ATAN(THH) 

IF(XB(K+1) .LT.XB(K)) TH( I ) =TH( I )+3 . 141592 

33 CONTINUE 

CALCULATE PANEL LENGTHS 

DO 34 I = N+l , M 
K = 1 + 1 

A( I ) = S0RT((XB(K+l)-XB(K))xx2+(YB(K+l)-YB(K))xx2) 

34 CONTINUE 

DO 35 I = N+l , M+l 

XB( I ) = XB( 1+1 ) 

YB ( I ) = YB( 1 + 1 ) 

35 CONTINUE 
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OOO ODD ono ooo oooo oooo 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

NON-LIFTING PART (2-BODIES) 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO 36 I = 1 , M 

CALCULATE ANGLE BETA AND NORMAL COMPONENT OF FREE STREAM 

VELOCITY — 



BEC I ) = TH(I)+3. 141592/2. 

VCI) = -(COS(BE(I))*COS(AN)+SIN(BE(I))xSINCAN)) 

— CALCULATE INFLUENCE COEFFICIENTS OF NORMAL VELOCITY — 

DO 37 J = 1,M 

IF(I.EQ.J) GO TO 37 

B = (XM(I)-XB(J))xCOS(BE(I))+CYM(I)-YB(J))xSIN(BE(I)) 

C = COS(TH(J))XCOSCBE(I))+SINCTHCJ))XSINCBECI)) 

E = 2.*C0S(TH(J))*(XM(I)-XBCJ))+2.xSINCTH(J))x(YM(I)-YBCJ)) 

F = (XM(I)-XB(J) )xx2+(YMCI)-YBCJ))xx2 
H = SQRT(4.XF-E**2) 

Z(I,J) = TEG(B,C,E,F,H,A(J)) 

7 CONTINUE 

ZCI,I)=3. 141592 
6 CONTINUE 

SOLVE SET OF LINEAR EQUATIONS FOR SOURCE STRENGTHS 

IDGT = 0. 

CALL LEQT2F (Z, 1 , M, 200 , V, I DGT , WKAREA, IER) 

— - CALCULATE I N FL U ENCE~COEFFI C IENT S OF TANGENTIAL VELOCITY 

DO 38 I = l.M 

DO 39 J = 1,M 

IF(I.EQ.J) GO TO 39 

B = (XM(I)-XBCJ))XCOS(TH(I))+(YM(I)-YBCJ))XSIN(TH(I)) 

C = COS(THCJ))xCOS(TH(I))+SIN(TH(J))XSINCTH(I)) 

E = 2.*C0SCTHC J))x(XM(I)-XB(J))+2.X3INCTH(J))X(YM(I)-YBCJ)) 
F = (XMCI)-XB(J))x*2+(YM(I)-YB(J))xx2 
H = SQRTC4.XF-EXX2) 

Y ( I , J ) = TEG(B,C,E,F,H,A(J)) 

39 CONTINUE 

Y(I, I ) = 0 . 0 

38 CONTINUE 

WRITEC8 , 95 ) 

CALCULATE TOTAL VELOCITY AND CP AT MIDPOINTS OF EACH PANEL 

DO 40 I = 1 , M 
S = 0 . 

DO 41 J = 1 , M 

S = S+V (J)XY(I,J) 

41 CONTINUE 

VV(I) = C0S(TH(I))XC0S(AN)+SIN(TH(I))XSIN(AN)+S 
CP(I) = 1.-VV(I)XX2 
WRITEC 8 , 93 ) I , XM( I ) , YM( I ) , VV( I ) , CP ( I ) 

40 CONTINUE 
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no o 



oooo noon oooo oooooooo 



ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

LIFTING PART (2-BODIES) 



CALCULATE SOURCE STRENGTHS DUE TO CIRCULATORY FLOW OF UNIT 

STRENGTH PAST FIRST BODY 

DO 42 I = 1 ,M 
S = 0. 

DO 43 J = 1,M 
G = 1. 

IFCJ.GT.N) G=0 . 

S = S+Y(I,J)*G 

43 CONTINUE 

V1(I)=S 

42 CONTINUE 
I DGT =0. 

CALL LEQT2F ( Z, 1 , M, 200 , VI , I DGT, WKAREA, IER) 

CALCULATE DISTURBANCE VELOCITY DUE TO CIRCULATORY FLOW OF UNIT 

STRENGTH PAST FIRST BODY 

DO 44 I = 1,M 
S = 0. 

DO 45 J = 1,M 
G = 1. 

IFCJ.GT.N) G=0 . 

S = S+YCI, J)*V1( J)+ZCI, J)*G 

45 CONTINUE 

VC1(I)=S 

44 CONTINUE 

--- CALCULATE SOURCE STRENGTHS DUE TO CIRCULATORY FLOW OF UNIT 

STRENGTH PAST SECOND BODY 

DO 46 I = 1/M 
S = 0. 

DO 47 J = 1,M 
G = 1. 

IF(J.LE.N) G=0. 

S = S+YCI, J)*G 

47 CONTINUE 

V 2 ( I ) = S 

46 CONTINUE 
I DGT=0 . 

CALL LEQT2F (Z,1,M,200,V2, I DGT , WKAREA, IER) 

CALCULATE DISTURBANCE VELOCITY DUE TO CIRCULATORY FLOW OF UNIT 

STRENGTH PAST SECOND BODY 

DO 48 I = 1,M 
S = 0. 

DO 49 J = 1,M 
G = 1. 

IFCJ.LE.N) G=0 . 

S = S+Y (I,J)XV2(J)+Z(I,J)XG 

49 CONTINUE 

VC2( I ) =S 

48 CONTINUE 
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non 



ono oor> non 



— CALCULATE VORTEX STRENGTHS BY KUTTA CONDITION — 

XI = VC1(1)+VC1(N) 

X2 = VC2C 1 )+VC2( N ) 

X3 = -VVd)-VV(N) 

XA = VCKN + 1)+VC1(M) 

X5 = VC2(N+1)+VC2(M) 

X6 = -VV(N+1)-VV(M) 

SI2 = (X3XXA-X1XX6 VCXAXX2-X1XX5) 

SI1 = CX3-X2*SI2)/X1 

CALCULATE TOTAL VELOCITY AND CP AT MIDPOINTS OF EACH PANEL — 

WRITEC8 , 96 ) 

DO 50 I = 1,M 

VT( I ) = VC1CI)XSI1+VC2(I)XSI2 
VVV(I) = VV( I )+VT ( I ) 

CP(I) = l.-VVV(I)**2 

WRITE(8,93)I,XMd),YMd),VVV(I),CPd) 

50 CONTINUE 
700 WRITEC6 , 97 ) 

97 FORMATdX, 'COMPUTATION COMPLETED') 

STOP 

END 

THIS FUNTION EVALUATES THE INTEGRALS (INFLUENCE COEFFICIENTS) 

FUNCTION TEG(B,C,P,Q,R,S) 

TERM1 = AL0GC(SXX2-P*S+Q)/Q) 

TERM2 = ATANC(2.XS-P)/R)-ATAN(-P/R) 

TEG = -CxTERMl/2.+(2.XB-Cxp)XTERM2/R 

RETURN 

END 
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IV. BOX METHOD 



A. INTRODUCTION 

The thin shear layer equations are more complicated than 
Laplace's equation because they are nonlinear. This chapter 
presents the box-method, which can be applied to the solution 
of the thin shear layer equations. The box method was intro- 
duced by Keller in 1970 [Ref. 4]. 

One of the basic ideas of the box method is to write the 
governing system of equations in the form of a first-order 
system. This system is solved by finite-difference approxima- 
tions and Newton's method is applied to solve the equations. 
Finally, the resulting linear system is solved by the block- 
elimination method. 

B. FALKNER-SKAN TRANSFORMATION 

The thin shear layer equations for incompressible laminar 
flow take the form 



ou , 3v 
8x dy 



0 



(4.1) 



u 





1 dP 
p dx 



+ v 



9 2 u 

9y 2 



(4.2) 



Boundary conditions are prescribed at the surface 



y = 0 u = 0 v = 0 
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and at the edge of the boundary layer 



y->°° u = (x) (4.3) 

It is convenient to reformulate the equations using the 
streamf unction and the similarity concept. Therefore the 
Falkner-Skan transformation is introduced. 

U 1/2 Re 1/2 

n = (— ) y = — - — v (4.4) 

VX 2 X 

ip(x,y) = (U q vx) 1//2 f (x, n) (4.5) 

Substituting Eqs. (4.4) and (4.5) into Eq. (4.2), we get the 
transformed momentum equation for 2-D laminar flows. 

f." + +m [i-(f») 2 ] = X (f||l-f-||) (4.6) 



where 



m 



U dx 
e 

(dimensionless pressure-gradient) 



with the boundary conditions 

n = 0 f ' = 0 f = 0 

n = n f ' = l 

1 1 oo 



(4.7) 
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If f is a function of n only, the right-hand terms of Eq. 

(4.6) will be zero. Then this will be a third-order ordinary 
differential equation whose solution is well-known as a 
"similar flow." But, if f is a function of n and x (non- 
similar flows), we need a numerical procedure, such as the 
box method. 

C. NUMERICAL FORMULATION (BOX METHOD) 

First of all, the coordinates (x,y) of a given geometry 
must be transformed to coordinates (£/h) to apply the Box 
method (see Figure 4.1). 



, Stagnati on Poi nt 




Figure 4.1. Transformed Coordinates of Upper Surface 
Airfoil and Net Rectangle for Difference 
Approximations 
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The boundary layer thickness in transformed coordinates is 
nearly independent of the streamwise distance and can be 
represented by a fixed number of profile points at fixed 
spacing . 

One of the basic ideas of the Box method is to write the 
governing system of equations in the form of a first order 
system. We write Eq. (4.6) in terms of a first-order system 
of PDE's 

f' = u(£,n) (4.8a) 

U ' = v (£ ,n ) (4.8b) 

(bv) ' + (5^p) fv + m(l-u 2 ) = £ (u - v |-|-) (4.8c) 

where a prime denotes differentiation with respect to n 
and b = 1 + v /v 

with the boundary condition 

f(£,0) = 0, u (£ , 0 ) = 0, u(£,nj = 1 (4.9) 

We denote the net points shown in Figure 4.1 as 

£q = 0 = £^_i + k^ i = 1 , 2 , . . . , I 

n 0 = 0 rij = r)j_^ + hj j = 1 , 2 , . . . , J — n ra 
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And we can introduce the following approximations: 



A) Coordinates of midpoints (£ ^,n and net functions 



1 * J j 

(g stands for f, u or v) 

? i = + S,-_i) n i e j(n- + n j -]L ) 



J 1^ 2 vs i T s i-l 

'2 



(4.10a) 



1-- D-o 



p * HS' 1 ’ vi 5 



1 / i . i 



(4.10b) 



where [ means the quantities (f or u or v) at point (£^,rij) 

B) Finite-difference approximation 

From Eqs . (4.8a) and (4.8b), the centered-dif f erence 

derivatives are 



f -f , 

_J Hi 

h . 

J 



= u 



. 1 

3 2 



(4.11a) 



l l 

u . - u . . 
J Hi 



= v 



. 1 

3 2 



(4.11b) 



After introducing these approximations into Eq. (4.8) and 

rearranging (the known quantities are moved to the right hand 

side), we get the equation (4.12c) which is centered about the 

point (5 i , n -^) . This represents the relationship of quan- 
i_ 2 3 ~2 

tities between the points of the box. 

h . 

d - - -f(uj + uj ,) = 0 (4.12a) 

j J - 1 - ^ J j 
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h . 



U i " U j-1 " ^ (v j " v j-l } " 0 



(4.12b) 



(b j V j - b j-l V j-l )/h j + a l (fv)i 1 - a 2 (u2)i 1 

3 2 3 ~ 2 



+ a (v -^f 

j ~2 3 ~2 



,i-l i 
" f . l v . 1 
3 2 3 2 



) = R 



i-1 
1 

3 ~1 



(4.12c) 



where 



a = 



' . 1 
1 2 



a. 



m . +1 

l 



+ ot , ot . 



m . + a 
i 



R 



i-1 
. 1 

3 ~2 



7 i-1 

{-L 1 +a[(fv) 1 ~(u ) 1 ]} -m 1 



3 ~2 



3 2 3 2 



i-1 
J . 1 
3 ~2 



b . v . -b . , v . , , , ~ i-1 

{ .J . 2—J 3 — + 5Ll(f v ) 1 +m[l-(u ) -J } 

3 3 ~2 



3 ~2 



The last of the above equations is non-linear. Therefore we 
introduce Newton's method to solve this system. We set 



(n+1 ) (n) . (n) 

g . + og . 

1 J J 



n = 0,1,2, 



(4.13) 



where the superscript in parentheses is the iteration counter 
with initial condition 
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( 0 ) 



= 0 ti 



( 0 ) 



= 0 v 



( 0 ) 



= v. 



i-1 



f : 

3 



( 0 ) 



-i-l (0) 
f . u . 



u i-i v (o> = v i-i 

3 3 3 



(4.14) 



(0) = -i-l (0) 

J J U J 



= 1 v 



( 0 ) 



= v. 



i-l 



1 < j < J-l 



Introducing Eq. (4.13) to Eq. (4.12) and dropping the quadratic 
terms in g^, we get (superscripts i and n are dropped for 
simplicity) 



Sf j ■ Sf j-i • 4 (Su j + 6u j-i> 


- 


(r i> j 


(4.15a) 


h . 

Su. - Su . , - — J-(Sv. + 5v . -i ) 
3 2 3 3-1 






(4.15b) 


(S,).6v. + (S,).6v . .+ (S^) - 5 f . + 


(s 4 ) 


. 6 f . , 
3 3-1 




+ (S c ) . Su . + (S,, ) . Su . , = 

5 ] 3 63 3 - 1 


( r 2 ) 


j 


(4.15c) 



where all terms are explained in Reference 5, with the 
boundary conditions 

Sf Q = 0 Su Q = 0 SUj = 0 (4.16) 

D. BLOCK ELIMINATION METHOD 

This is a very effective way to solve linear difference 
equations, discussed by Keller in 1974 [Ref. 6]. We write 
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Eq. (4.15)- in a matrix-vector form 



A6 = 



(4.17) 



where 



A = 



A o c o 




f \ 

4 o 




r o 


B. A. C. 










Ill 




1 




r.. 


X 








1 


s'* 

\ s v 


6 2 


• 


r = 


• 






• 




• 


v ' 










B t A -0-i 










j-i j-i j-i 










A- 




6 






J J 




J 




r J 



where 



A„ = 



1 


0 


0 




1 


-V 2 


0 


0 


1 


0 


A J E 


(S 3>J 


(S 5»J 


(S 1>J 


0 

l. 


-1 


-hj/2 




0 

L 


1 


0 





1 -h./2 

3 


0 




0 


0 


0 


A. = 

3 


(s 3 ). (s 5 ). 


(s l»j 


c. = 

3 


0 


0 


0 




0 -1 


-Vl 




0 


1 


-v / 2 



1 < j < J-l 



B. 

D 



-1 -h./2 

3 

(S 4 )j (S 6)j 



«2>j 



1 < j < J 



where 
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6 . 
D 



0 < j < J 



<5f . 

J 

6 u . 

J 

6 v . 

3J 



r 

J 



<r i > j 

(r 2>j 
(r 3> j 



1 < j < J-l 



(r l>J 
(r 2* J 



0 

0 



According to Keller's block elimination method, we have to 
factorize the matrix A. 



A = P xq 



(4.18) 



where 
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I is the identity matrix 







1 


0 0 










I 


— 


0 


1 0 














0 


0 1 










From Eq . ( 4 . 18 ), we find 


that 










o 

a 


= A o 










( 4 . 19 a) 


p j Q :-i 


= B . 
3 




j 


= 12 

J- r ^ t • • 




( 4 . 19 b) 




= A . - 
3 


Vj-i 


j 


= 12 

^ r ^ t • • 


• ,J 


( 4 . 19 c) 


Keller showed that 


the 


matrix P . 


has 


the same structure 


as 








3 










the matrix . From Eq. 


( 4 . 19 ) , 


we 


derive the elements 


of 


P . and Q . . 
















3 3 


















(Pn)j 


( ?12 } 


j 


(p 13 > j 






P . = 

3 


( P21>j 


( P 2 2 } 


j 


( P 2 3 > j 








0 




0 




0 

J 








(q ll } j 


( q i2 ) 


i . 

3 


j 






Q j 


(q 21 ) j 


^ q 22 ^ 


i . 

3 




0 < j < 


J-l 




0 




-1 




-*Vi /2 
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Each element of and is explained in Reference 5. 
we let 



Q& = W 



Introducing Eqs . (4.20) and (4.18) into Eq. (4.17), 



PW = r 



Then, from Eq. (4.21), we find that 



w o = r 0 



W . = r.-P.W.-, l<j<J 

3 3 ID" 1 - J “ 



The elements of are listed in Reference 5. Finally, 
get the matrix form to get 6 from Eq. (4.20) . 



lO 

o 

O 

o 




/ 

o 

V . 




V 


N 




1 — 1 

*0 




W 1 


x V 




• 




• 


\ 




I 




1 


1 — 1 

1 

l^> 

u 

1 — 1 

1 

o 










Q J 








W J 


■*«. i 




S J 







From Eq. (4.22), we find that 



Q J 6 J = W J 



If 



(4.20) 



(4.21) 



we 



(4.22) 



(4.23a) 



77 



Q j 6 j = W j " c j 6 j +i • j = . . . ,0 (4.23b) 



Therefore 6^ can be obtained by calculating the terms Q. ^ 

and W. (see Reference 5). 

3 
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V. INTERACTION METHOD 



A. INTRODUCTION 

Interactive methods provide a special coupling between 
viscous and inviscid flows. They are intended to compute 
flows including separation. Thus these methods may be regarded 
as an alternative to the Navier-Stokes solvers, which are 
hardly engineering tools because of their huge computer time 
and storage requirements. 

The classical method to compute viscous flows past airfoils 
proceeds as follows: The velocity distribution is computed 

by any appropriate inviscid flow solver. The output of the 
inviscid flow solver becomes the input of the viscous flow 
solver. Solving for viscous flow consists of integrating the 
boundary layer equations. Provided that the flow remains 
attached this procedure allows a reliable prediction of lift 
and drag. Information is transferred only once from inviscid 
to viscous regions. However, many flows cannot be modelled by 
one-time information transfers. 

Separation bubbles and separated flows especially require 
a close coupling between viscous and inviscid regions. Inter- 
action schemes provide a better exchange of information between 
these two regions. 

The elements of interaction schemes are: direct or inverse 

inviscid flow solver and direct or inverse viscous flow solver. 
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•Boundary Condition 
Flow 


Direct 


Inverse 


Inviscid 


• Zero normal 
velocity at 
the surface 


• Prescription of 
velocity 
distribution 


Viscous 


• No slip 
condition 

• Prescription of 
external 
velocity 


• No slip condition 

• Prescription of 
displacement 
thickness 



The direct boundary layer method has the disadvantage that 
the boundary layer equations are singular at the point of 
separation. However, if the external velocity is computed 
by prescribing a displacement thickness (known as the inverse 
boundary layer method) , they can be integrated through the 
point of separation. 

The next problem associated with the regions of reversed 
flow is numerical instability, because downstream marching 
procedures cannot be applied in regions of reversed flow. 

The most common approximation to get this instability under 
control, the so-called FLARE approximation, neglects the 
momentum transport term u 3u/9x in regions of reversed flow as 
long as this region is small. However, as the size of this 
region increases, the FLARE approximation becomes inaccurate. 

One of the procedures which can be taken into account is 
called the DUIT (Downstream, Upstream iteration) . It consists 
of a sequence of alternating up and downstream sweeps. 

There are several types of recently developed interaction 
models. All procedures have to solve both the inviscid (Laplace 
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equation) and viscous flow, whose equations can be written 
according to 



3u 3_v 
3x 3y 



0 



(5.1) 





(5.2) 



where 



b = l = constant in laminar flow 



b - 1 + v t /v in turbulent flow 



Four interaction models can be distinguished: Direct, Inverse, 

Semi-inverse, and Simultaneous interaction methods which are 
subject to different boundary conditions. 

The most advanced interaction scheme is the simultaneous 
interaction. We call it the "strong interaction" (direct and 
inverse interactions guarantee weak coupling only) . Examples 
in Section V.D are computed by the Cebeci program using this 
method. Good agreement is obtained between the results of 
interaction schemes and experimental results . 

B. FOUNDATION OF THE INTERACTION SCHEMES 

1. Direct Interaction Scheme 



inviscid and a direct viscous flow solver (see Figure 5.1a). 
The usual sequence is: 



The direct interaction model is composed of a direct 
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a) DIRECT 



Airfoil 

geometry 




c) SEMI- INVERSE 




b) INVERSE 



Airfoil 

geometry 

r 

INVISCID 
(Direct ) 



i Initial guess of 
j ! displacement thick 
v |.ness distribution 




Figure 5.1. Global Organization of Interaction Methods 
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(1) Calculate the external velocity distribution by 
inviscid flow computations. 

(2) Calculate the displacement thickness, 5*, by viscous 
flow computations using the external velocity as 
boundary condition. 

(3) Compute an updated shape of the displacement body and 
repeat steps 1 and 2 until the results converge. 

Unfortunately, the direct boundary layer method suffers 

numerical breakdown at the point of separation (Goldstein 

singularity) . Therefore this scheme is not appropriate to 

handle airfoil flows with separation. 

2 . Inverse Interaction Scheme 

This method was introduced to overcome the singularity 
problems near separation. It combines an inverse inviscid and 
an inverse viscous flow solver (see Figure 5.1b). However, 
the overall procedure suffers from very slow convergence. 

Due to this shortcoming one shall apply this inverse scheme 
to regions of separated flow only. 

3 . Semi-Inverse Interaction Scheme 

This method combines a direct inviscid flow solver 
with an inverse viscous flow solver with the same input (dis- 
placement thickness) . This leads to two external velocity 
distributions, U e ^(x) and U gV (x) (see Figure 5.1c). Satis- 
factory convergence is ensured by a relaxation formula, which 
is introduced to define an updated displacement thickness 
distribution . 



6 * 

new 



(x) 



U ev (x > 

4 Jld' x)I1 + -'57(51 

el 



D] 



(5.3) 
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where w is a relaxation parameter. The numerical weakness of 
the purely inverse scheme is improved by this method, but 
both inviscid and viscous regions are still coupled loosely. 

4 . Simultaneous Interaction Scheme 

The simultaneous interaction scheme emphasizes strong 
interaction between the outer inviscid and the inner viscous 
region. The external velocity U^tx) and the displacement 
thickness 6*(x) are treated as unknown quantities. An addi- 
tional relation is added, the so-called interaction law which 
can be given by the "blowing velocity" concept. 

The equations are solved by the inverse method with 
successive sweeps over the airfoil surface (see Figure 5. Id). 
This method is compatible with the weak interaction scheme 
where both inviscid and viscous regions are coupled loosely. 

For each sweep, the external velocity for the boundary 
layer equation is written as 

U(x) = U°(x)+6U(x) (5.4) 

e e e 

where 

U°(x) is the inviscid velocity; 

6U (x) is the perturbation due to the dis- 
placement effect of a boundary layer. 

The blowing velocity concept is introduced to get the pertur- 
bation velocity SU e by the interaction law. The displacement 
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effect, of a boundary layer can be modelled by ejecting fluid 
at the airfoil's surface (see Figure 5.2). 




Figure 5.2. Blowing Velocity Concept 



A properly arranged source distribution on the surface dis- 
places the streamlines away from the surface such that the 
virtual displacement body becomes a streamline. 

Our first, goal is to determine the source strengths such 
that the tangential flow condition on the displacement body 
takes the form 



v (x , 6 *) _ d6 * 

U (x) dx 

e 



To achieve this goal we use the thin airfoil approximation: 

(1) The displacement thickness is assumed to be so 
small that u-velocity components do not vary across 
the layer. 

(2) The airfoil in this connection can be represented by 
a straight line. This approximation implies that 
the blowing velocity v (x,0) equals half of the 
source strength. 
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Therefore , 



o (x) 
2 



v (x , 0) 



6 ^ 



3 v 



= v(x,s*) - / 



= u 



d6 



e dx 



+ 



dU 

e 

dx 



6 * 






(5.6) 



where ^(U e <5*) is defined as blowing velocity. Our second 
goal is to relate the perturbation velocity, <5U e to the blowing 
velocity. This process is quite similar to evaluating tangen- 
tial velocities in the panel method. In fact, this is even 
simpler because of the straight line surface. 



6 U 



e 




(Cl 

- C 



(5.7) 



where the interaction region is limited to a finite range 
x a _< x _< x^ . This integral is referred to as the Hilbert 
integral. Rewriting Eq . (5.4)/ we finally obtain the inter- 

action law. 



U (x) 
e 



x. 



U°(x) 



+ 






6 *) 



dC 

x -C 



(5.8) 
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The numerical implementation of Eq . (5.8) requires a discrete 
approximation of the Hilbert integral. This can be performed 
by using the trapezoidal rule. 



The examples in Section V.D demonstrate that this inter- 
action method can give reliable results for flows up to high 
angle of attack, including flows with bubbles and separation. 

C. CONSIDERATION OF BOUNDARY LAYER TRANSITION AND OF 
TURBULENT FLOW MODELLING 

1. Transition 



drag and lift of an airfoil is the transition point. Boundary 
layer transition is affected by many parameters, for example, 
the pressure distribution (major parameter) , the wall roughness 
and the intensity of the free stream turbulence, etc. Because 
of this fact, the theoretical modeling of transition is very 
complicated and one therefore resorts to experimental 
information . 



correlation formula is used, which was given by Cebeci and 
Smith (1974) as a relation between R and Re at transition. 



One of the most important parameters to predict the 



In the Cebeci program, the following experimental 



x 



R, 



6 



tr 




(5.9) 



where 



Re 



x 



U x/v 
e 
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and 0 is the momentum thickness. 

2 . Turbulent Flow Model 

Unlike laminar flows, turbulent flows have a compli- 
cated time-dependent behavior. It is too difficult to deal 
with the instantaneous properties. Thus, the mean-flow 
properties are applied in turbulent flow. 

The most common mean-flow models are the "eddy- 
viscosity" formula which are based on thin shear layer 
assumptions . 



where is related empirically to the mean flow velocity 
gradient and the length scale. In the Cebeci program, v 
is presented by the algebraic eddy-viscosity formulation of 
Cebeci and Smith. 



pu V 




(5.10) 



! 



(0.4y [1-exp (-y/A) ] } 2 ||^|y tr 0 £ y £ y c 



V 



t 



(5.11) 



oo 



“I / (u -u)dy|y y 
0 



y < y < 6 

J c — J — 



More detailed descriptions are listed in Reference 7. 



D. EXAMPLES 



The subsequent examples were computed using a program 
developed by Cebeci and coworkers [Ref. 17] , on the NPS IBM 37 0. 



1 . 



Demonstration of the Program Capabilities 



The velocity profiles on both upper and lower surfaces, 
as well as in the wake, are presented in Figure 5.3. At 
this angle of attack (a = 10°) , transition occurs very close 
to the leading edge on the upper surface. The boundary layer 
thickness is quite thin in the accelerated flow region (right 
after the leading edge) , but it grows thicker farther down- 
stream in the decelerated flow region (near the trailing edge) . 
Eventually, we find a small reversed flow region just before 
the trailing edge in this case. The wake region shows the 
mixing layer which decays with increasing downstream distance. 

Figure 5.4 demonstrates how lift, drag and the loca- 
tion of transition depend on the angle of attack. The skin 
friction drag is dominant at low angles of .attack, the pressure 
drag at high angles of attack (see Figure 5.4b). 

Figure 5.5 shows the distributions of the skin friction 
coefficient, displacement thickness and shape parameter in 
dependence of Reynolds number and angle of attack. In the 
attached flow, the skin friction coefficient decreases along 
the downstream direction until the point of transition 
(laminar region), but increases steeply after transition and 
then decreases again because the skin friction is related to 
the slope of the velocity profile, 3u/3y, at the surface. At 
high angles of attack, transition on the upper surface occurs 
close to the leading edge and the negative skin friction 
coefficient near the trailing edge indicates separated flow. 
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Fi gure 5.3. 



Velocity Profiles on 
of FX 63-137 Airfoil 



the Surface and Wake 
(a = 10° , Re = 5 > 10 5 ) 
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Figure 5.4, (a) Lift, (b) Drag, (c) Transit 

Calculations (FX 63-137 Airfoil, 
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(b) Constant angle of attack , 2 Reynolds Number ^ attached flow(5*10 6 ) • -- 

( 3° ) a short bubble ( . Cx 10 s ) : 

Figure 5.5. Skin Friction Coefficient, Displacement 

Thickness and Shape Parameter for FX 63-137 



Also the displacement thickness and the shapeparameter increase 
drastically in regions of separated flow (see Figure 5.5a). 

At low Reynolds numbers, laminar flow can separate at 
mid-chord and reattach (we call this phenomenon a "bubble"). 
Reattachment can often occur if the pressure gradient de- 
creases rapidly soon after separation, so that a strong reversed 
flow is not established. Thus the shear layer reattaches onto 
the surface. Accordingly, the displacement thickness and the 
shape parameter increase in the bubble region (see Figure 
5.5b) . 

2 . Comparison with Experimental Results 

The comparison of pressure distributions is shown in 
Figure 5.6. The overall lift of inviscid calculations devi- 
ates 20% from the experimental results. The interaction method 
improves the accuracy, but the computation still overestimates 
the lift. The lift coefficient obtained by Cebeci’s program 
is approximately 10% larger than the measured one (see Reference 
8 ) . 

Figure 5.7 demonstrates that the accuracy of this method 
drops with lower Reynolds numbers. The reasons for this de- 
creased accuracy at low Reynolds numbers are the higher proba- 
bility to separate and the used turbulence model, which seems 
to be inappropriate for low Reynolds numbers. Therefore, low 
Reynolds numbers and high angles of attack pose severe limi- 
tations (see Figure 5.7a & b) . The experimental results are 
taken from References 7 and 9. 
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Figure 5.6. Comparison of Pressure Distributions (HACA 

4412 Airfoil, a = 12.15° 
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; Experimental 

: Interaction method 

Figure 5.7 C 9 and C d Curves for (a) FX 63-137, Re = 199487 
(b) GA-1 , Re = 1900000, (c) FX 63-137, Re = 280000 
( Source of experimental results: (a)Ref.9 (b)Ref.7 (c)Ref. 10 



On the other hand, experimental measurements should 



not be expected to be exact. Turbulence level and interference 
effects influence the wind tunnel measurements. Figure 5.7c 
shows good agreement between computed and experimental results 
taken from Reference 10, at low Reynolds numbers. 

The location of transition and separation points have 
an important influence on the lift and drag coefficients. 

Figure 5.8 shows that the prediction of laminar separation, 
reattachment, transition and separation points on the airfoil 
surfaces are in reasonable agreement with the experimental 
data taken from Reference 9. 



96 



VO 

O 




- w 



UJ 

© 

a 

a 

< 

HI 

< 



CO 







4-* 


o 


OJ 

Cl 


-C 


— " 


g 






ro 


c 


4-> 


o 


C 


•f” 


OJ 


«LJ 


cE 


o 


•f- 


t tJ 


S- 


s- 


<u 


o; 


Cl 




X 


c 


LlJ 


*— « 



X 



o/x 



c 

o 




w 

U 3 



a 

u? 

c 

< 



-2 

< 



97 



Figure 5.8. Transition and Separation Positions for NACA 64 3 -418 



VI. SUMMARY AND RECOMMENDATIONS 



This thesis treats the problem of incompressible two- 
dimensional steady flow past airfoils or airfoil combinations 
at large angles of attack. A panel method was developed to 
compute the inviscid flow over two cylinders, airfoil-flap 
combinations and airfoils in ground effect. In addition, 
Cebeci's viscous/inviscid interaction method was applied to 
several airfoils and compared iwth available experimental 
data. The agreement is generally quite encouraging. The 
calculations show the sensitivity of the computations to 
Reynolds number and transition. More work is required to 
evaluate the potential and limitations of the viscous/inviscid 
interaction method. 
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